Bifurcations and stability of gap solitons in periodic potentials 
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We analyze the existence, stability, and internal modes of gap solitons in nonlinear periodic 
systems described by the nonlinear Schrodinger equation with a sinusoidal potential, such as pho- 
tonic crystals, waveguide arrays, optically-induced photonic lattices, and Bose-Einstein condensates 
loaded onto an optical lattice. We study bifurcations of gap solitons from the band edges of the 
Floquet-Bloch spectrum, and show that gap solitons can appear near all lower or upper band edges 
of the spectrum, for focusing or defocusing nonlinearity, respectively. We show that, in general, two 
types of gap solitons can bifurcate from each band edge, and one of those two is always unstable. A 
gap soliton corresponding to a given band edge is shown to possess a number of internal modes that 
bifurcate from all band edges of the same polarity. We demonstrate that stability of gap solitons 
is determined by location of the internal modes with respect to the spectral bands of the inverted 
spectrum and, when they overlap, complex eigenvalues give rise to oscillatory instabilities of gap 
solitons. 

PACS numbers: 42.65.Tg, 42.65.Jx, 42.70.Qs, 03.75.Lm 



I. INTRODUCTION 

Periodic structures are very common in physical prob- 
lems, with the crystalline lattice being the most famil- 
iar classical example. One of the important features of 
such systems is the existence of multiple frequency gaps 
in the wave transmission spectra. Such spectral gaps arc 
responsible for a strong modification of the wave disper- 
sion and diffraction that occurs when waves experience 
resonant Bragg scattering from a periodic structure [1]. 

When nonlinear self-action becomes important, the 
systems with periodically modulated parameters demon- 
strate a number of new effects; in particular they can sup- 
port a novel type of solitons, the so-called gap solitons, 
which exist in the gaps of the linear wave spectrum due to 
a strong Bragg scattering and coupling between the for- 
ward and backward propagating modes [2, 3, 4]. During 
the last years, it was shown that gap solitons may exist 
in different types of nonlinear periodic structures includ- 
ing low-dimensional photonic crystals and photonic lay- 
ered structures [5, 6], waveguide arrays [7, 8], optically- 
induced photonic lattices [9, 10], and Bose-Einstein con- 
densates loaded onto an optical lattice [11, 12, 13]. 

There are known two standard approaches to study 
nonlinear localized modes and gap solitons in periodic 
structures [14]. The first approach is based on the deriva- 
tion of an effective discrete nonlinear Schrodinger equa- 
tion and the analysis of its stationary localized solutions 
in the form of discrete localized modes or discrete soli- 
tons [15]. In the solid-state physics, the similar approach 
is known as the tight-binding approximation that, in ap- 
plication to periodic photonic structures, corresponds to 
the case of weakly coupled defect modes excited in each 
individual waveguide of the structure. The analogous 



concepts appear in the study of other systems such as 
the Bose-Einstein condensates in optical lattices [16]. 

On the other hand, weak nonlinear effects in optical 
fibers with a periodic modulation of the refractive in- 
dex are well studied in the framework of the other famil- 
iar and well-accepted approach, the coupled-mode the- 
ory [4] . The coupled-mode theory for nonlinear periodic 
structures is based on a decomposition of the wave field 
into the forward and backward propagating modes, under 
the condition of the Bragg resonance, and the deriva- 
tion a system of coupled nonlinear equations for those 
two modes. Such an approach is usually applied to an- 
alyze nonlinear localized waves in the systems with a 
weakly modulated optical refractive index known as gap 
(or Bragg) solitons. 

A number of recent experiments in the nonlinear 
guided wave optics [7, 8, 10, 17, 18] and Bose-Einstein 
condensates [13, 19] were conducted in periodic struc- 
tures under the conditions when none of those approx- 
imations is valid. Indeed, one of the main features of 
the wave propagation in periodic structures is the exis- 
tence of a set of multiple forbidden gaps in the trans- 
mission spectrum. As a result, the nonlinear ly-induccd 
localization of waves can become possible in each of these 
gaps [4, 7, 20]. However, the effective discrete equation 
derived in the tight-binding approximation describes only 
one transmission band surrounded by two semi-infinite 
band gaps and, therefore, a fine structure of the band- 
gap spectrum associated with the wave transmission in 
periodic media is lost in this approach. On the other 
hand, the coupled-mode theory of gap solitons describes 
only an isolated narrow gap in between two semi-infinite 
transmission bands, and it does not allow to consider si- 
multaneously the localized modes due to the total inter- 
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nal reflection as well as to study the band coupling and 
inter-band resonances. Recently, it was realized that the 
study of the simultaneous existence of localized modes of 
different types is a very important issue in the analysis 
of stability of nonlinear localized modes and gap soli- 
tons [20, 21]. 

The main purpose of this paper is to develop, for the 
first time to our knowledge, a general analytical descrip- 
tion of the existence, bifurcations, and stability of spa- 
tially localized nonlinear modes (i.e., lattice and gap soli- 
tons) in the framework of an effective continuous model 
described by the nonlinear Schrodinger equation with a 
periodic external potential. The use of this well-accepted 
nonlinear model for our analysis allows us to remove all 
restrictions of both approaches mentioned above, and to 
study consistently the effects of the bandgap spectrum 
on the properties and linear stability of gap solitons. 

First, by applying the multi-scale asymptotic analyti- 
cal methods, we show that such gap solitons may appear 
in all band gaps of the periodic potential for any sign of 
nonlincarity, but they bifurcate from different band edges 
for different signs of nonlincarity. Second, we demon- 
strate that, in general, only two branches of gap solitons 
bifurcate from each band edge, and one of those two is 
always linearly unstable. Third, we study stability of gap 
solitons in a selected band gap and find the soliton inter- 
nal modes[22, 23] bifurcating from all other band edges of 
the same polarity. However, only one internal mode can 
bifurcate from the band edge where the gap soliton orig- 
inates itself. At last, we analyze the conditions when the 
bifurcation of the internal modes can give rise to complex 
eigenvalues, which are shown to be responsible for oscil- 
latory instabilities of gap solitons due to single-gap [24] 
or inter-gap [20, 21] resonances. 

The paper is organized as follows. Section II presents 
our physical model which is described by an effective non- 
linear Schrodinger equation with an external periodic po- 
tential of the simplest sinusoidal form. Section III sum- 
marizes the studies of the spectral properties of the linear 
eigenvalue problem with a periodic potential. In Sec- 
tion IV we study bifurcations of gap solitons by means of 
the weakly nonlinear approximation. Section V presents 
the analysis of the exponentially small corrections beyond 
the weakly nonlinear approximation. Section VI dis- 
cusses the stability problem of gap solitons. Symmetry- 
breaking instabilities are studied in Section VII. Inter- 
nal and oscillatory instability modes of gap solitons as- 
sociated with non-zero eigenvalues are studied in Sec- 
tion VIII. Finally, Sec. IX summarizes our results and 
discuss further perspectives. Appendix A gives details 
of the numerical method for calculations of eigenvalues. 
Appendices B and C present details of derivations, which 
are used in Sections VII and VIII, respectively. 



II. MODEL 

We consider the cubic nonlinear Schrodinger (NLS) 
equation with an external periodic potential in the form, 

m = -*xx + V(x)* + a|*| 2 f, (1) 

where V{x + d) = V{x), d is the fundamental period, 
and a = ±1 defines the type of the wave self-action 
effect, namely self-focusing (er = —1) or self-defocusing 
(a = +1). The analytical results presented below are 
rather general, and they are valid for different types of 
smooth arbitrary-shaped periodic potentials. However, 
in the numerical examples discussed below we consider 
the squared sine potential, 

V(x)=V sm 2 {^). (2) 

The harmonic potential (2) describes, in the mean-field 
approximation, the dynamics of the Bosc-Einstein con- 
densate in an optical lattice, when the parabolic trap is 
neglected [11, 12]. The squared sine potential V(x) has 
two extremum points on the period of x, such that x = 
is the minimum and x = d/2 is the maximum of V(x). 

Stationary localized solutions of the cubic NLS equa- 
tion (1) for gap solitons are found in the form ^S(x,t) = 
ip(x) exp(— i/ii), where \x is referred to as the soliton pa- 
rameter. The soliton profile ip(x) is found as a spatially 
localized solution of the nonlinear problem: 

-ip" + V(x)^ + a\^\ 2 ip = flip, (3) 

where the prime stands for a derivative in x. Existence 
and multiplicity of multi-humped localized states ip{ x ) m 
the spectral gaps of the periodic potential V{x) were con- 
sidered by means of the variational methods by Alama 
and Li [25, 26]. Bifurcations of bound states were ana- 
lyzed by Kupper and Stuart [27, 28] and Heinz and Stu- 
art [29] who proved that, depending on the sign of the 
nonlinear term, lower or upper end-points of the contin- 
uous spectrum are bifurcation points. Extension to the 
multi-dimensional case was developed with Bloch waves 
of the linear Schrodinger operator [30]. The number of 
branches of bound states was classified in terms of the 
eigenvalues of the linear Schrodinger operator with pe- 
riodic and decaying potentials [31]. Eigenvalues of the 
latter (linear) problem were previously considered by 
Gesztesy et al. [32] and Alama et al. [33]. In applica- 
tion to the problem of the Bose-Einstein condensates in 
optical lattices, the stationary model (3) has been consid- 
ered recently by Louis et al. [11] who found numerically 
different types of spatially localized solutions in different 
band gaps. 

All previous results were restricted to the study of the 
existence of spatially localized solutions. Here, we use 
more general methods (but, in some sense, less rigorous 
from the mathematical point of view) and study bifur- 
cations, stability, and internal modes of gap solitons. To 
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FIG. 1: The structure of spectral bands of the linear periodic 
problem (4): (a) Trace of fundamental matrix A versus /i; (b) 
Floquet exponent k versus fi; (c) Solid: Bloch waves at the 
band edges, as indicated by arrows; Dashed: potential V(x). 
Parameters are Vb = 1 and d = 10. 

achieve these objectives, we apply the multi-scale per- 
turbation series expansion methods, developed earlier by 
Iizuka [34] and Iizuka and Wadati [35]. With the per- 
turbation series methods, we classify systematically the 
branches of gap solitons bifurcating from the band edges 
to the band gaps, as well as their stability. 

III. SPECTRAL BANDS AND GAPS 

Periodic potential V(x) induces a band-gap structure 
in the linear Schrodinger spectral problem: 

-V>" + V(x)ij) = flip. (4) 

The spectral bands are located for fi G Sband, where we 
enumerate the band edges in the following order: 

Sband = [/"0,A*l] U [/i3,M2] U [fin, fis] U [fi 7l fie]... (5) 

The spectral bands are computed for the squared sine po- 
tential (2) and the results are shown on Figure 1. Review 
of spectral theory for periodic potentials can be found in 
the book by Eastham [36] . Here we recover some details 
which arc important for our analysis. 

When the spectral parameter fi is taken inside the 
spectral bands, i.e. fi £ Sband, the problem (4) has two 
linearly independent solutions in the form of Bloch waves, 

Vi = M^> ik(fl)x , ^ - Mx)e~ ikWx , (6) 

where 4>i,2{x) are periodic functions and k(fi) is the Flo- 
quet exponent, which can be chosen inside the first Bril- 
louin zone such that < k(fi) < ^. The graph of k(fi) is 
shown in Fig. 1(b) for the first three spectral bands. 

The spectral bands of the periodic potential (4) are 
described by the function A(/i), which is the trace of 
fundamental matrix of solutions [36] 

A(fi) = 2cosk(fi)d, (7) 



The spectral bands are defined for —2 < A(/i) < 2, which 
corresponds to propagating waves with real k. On the 
other hand, the waves become exponentially localized in- 
side the gaps, where |A(//)| > 2 and Im(fc) 7^ 0. A char- 
acteristic dependence A(fi) is displayed in Fig. 1(a). 

There are infinitely many spectral bands for a one- 
dimensional periodic potential V(x), where |A(/i)| < 
2 [36]. If A' (fi n ) 7^ at the band edge fi = fi„, two 
adjacent spectral bands do not overlap, such that the cor- 
responding band gap has a non-zero width. We consider 
the non-degenerate spectral band, such that A'(fi n ) 7^ 
at the end point fi = fi n . 

The even- numbered band edges fi — fi2mi rn > corre- 
spond to periodic Bloch functions ip2m(x + d) = ip2m(x), 
while the odd-numbered band edges fi = fi2m—it m > 1 
correspond to anti-periodic Bloch functions ip 2 m-i(x + 
d) = —i))2m—i(x)- The Bloch functions ip n (x) for the 
first five band edges fi = fiQ, fix, 113, fi2, fii are shown on 
Fig. 1(c). 

We now demonstrate that the bifurcations of bound 
states and stationary gap solitons may occur when the 
two fundamental solutions ipi,2(x) in (6) become linearly 
dependent. Since 4>i(x)e 2lk ^ x solves the same equation 
as 4>2(x) but it is not a periodic function of x, unless 
k(fi) = (mod( 2 f)) or k(fi) = § (mod(^)), the two 
solutions ^1,2(2:) are always linearly independent in the 
interior of the spectral bands fi <E Sband- On the other 
hand, the two solutions ipizfa) become linearly depen- 
dent at the band edges fi — fi n , since 

1p2m(x) = <t>l(x) = <h{n) ( 8 ) 

and k(fi2m) = (mod (^f )) at the even- numbered band 
edges, and 

4> 2m -i(x) = (/)i(x)e 2L ^' = (j) 2 (x)e~z~ (9) 

and k(fi2m-i) = 2 ( moc K^f)) at tne odd-numbered 
band edges. The band edge fi = fi n has geometric multi- 
plicity one with the only linearly independent Bloch func- 
tion ip n (x). The second, linearly independent solution of 
(4) at fi = fi n grows linearly in x. The band edge fi = fi n 
has however algebraic multiplicity two, since there ex- 
ists a generalized Bloch function ip„ (x) that solves the 
derivative problem: 

+ V(x)^ - fi n ^ = 2^' n (x). (10) 

It follows from (10) that the generalized Bloch functions 

^2m( x ) an d i , 2m-\( x ) are P er i°dic and anti-periodic in 
x, respectively. We conclude that the band edges fi = fi n 
are the only bifurcation points of the linear spectrum 
M £ Sband, associated with the periodic potential V(x). 

The band curvature near the band edge fi — fi n follows 
from the expansion of A(fi) defined in (7): 

A(fl n ) + A'(fl n )(fl - fin) + 0(fl - fl n ) 2 

= (-1)" (2 - d 2 (k - k n ) 2 + 0(k - k n ) 4 ) , (11) 
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where fem = and fem-i = 5 - As a result, A(/x n ) 
2(-l) n and 



m = fi n - $Xk - k n f + o(fc - k n y 



(12) 



such that 



,(2) 



d 2 (-l) n 
A'(/i„) 



(2) 

The band curvatures can be expressed in terms of 
Bloch functions ip n (x) and ipn\x) at /i = [i n . We use 
perturbation series expansions for fundamental solutions 
4>\^{x) near the band edges \x = fi n : 

d>i,2{x)e ±ik ~ x = ij n (x) ±i(k- k n )^\x) 

-(k-k n f^\x) + 0{k-k n f. (13) 

The eigenvalue [i is expanded in the perturbation series 

( 2) 

(12). The second-order correction ipn (x) satisfies the 
non-homogeneous linear equation: 

-{^)" + V{x)^-^ 

= 2(v4 1) )'+(l + Mi 2) )V' ra - (14) 

If 4>\.2{x) are periodic functions of x, the second-order 

correction tpn \x) in the perturbation series (13) is peri- 
odic for n = 2m and anti-periodic for n = 2m — 1. By 
Frcdholm Alternative, this condition is satisfied if the 
right-hand-side of (14) satisfies the constraint: 

(l + M £ 2) ) J i'ldx + 2 J V„ tyP) ' dx = 0. (15) 

Therefore, the band curvature /A is expressed in terms 
of integrals of ip n (x) and tpn\x). The perturbation sc- 
ries expansions (12) and (13) can be continued algorith- 
mically to the higher orders in powers of (k — k n ). 



IV. BIFURCATIONS OF GAP SOLITONS 

Nonlinear bound states (gap solitons) of the NLS equa- 
tion (1) are stationary solutions in the form: 



(16) 



where the real function $ s (x) decays to zero as |x| — > oo 
and satisfies the nonlinear problem, 



(17) 



When Q s {x) is small, the nonlinear potential a§ 2 s (x) acts 
as a perturbation term to the periodic potential V(x). 
The perturbation term leads to bifurcation of gap solitons 
<f> s (x) from the band edges /i = fi n of the linear band-gap 



spectrum. We study bifurcations of gap solitons with the 
multi-scale perturbation series expansions: 



H, = /i„ + e 2 A„ + 0(e 4 ), 



(18) 



and 



$ s (x) = e<f> t (x;X), X = e(x-x ), e < 1, (19) 
where 

$ e (x\X) = A n {X)^ n {x)+eA' n {X)^\x) 

+ e 2 4>?Hx;X) + 0(e 3 ). (20) 

Here A n {X) is a space-decaying bound state and ipn(x) is 
the periodic or anti-periodic Bloch function. Parameter 
xq determines a location of A„(X) with respect to ij) n (x). 
The Bloch functions ip n {x) and ip£\x) are defined from 
the linear problems (4) and (10). 

(2) 

The second-order correction term (f> s (x; X) satisfies 
the linear non-homogeneous equation: 

-^)y+v(x)^-^ s 2) 

= A'; U n + 2 (Vi 1 ^) - + A«A„^, (21) 

The secular growth of <fi s (x; X) in x is removed if the 
right-hand-side of (21) satisfies the Frcdholm condition, 
which reduces to the nonlinear equation for A n = A n (X): 



where 



u (2) All , (2) .3 _ A A _q 



(2) _ Jo fPtdx 

An u r d 



Jo ^Idx 



(22) 



(23) 



Using the constraint (22), we represent the second-order 

( 2) 

correction term <f) s (x;X) in the form: 

^Hx;X)=A:(X)^\x)+AI(X)^ 12 \x), (24) 

(2) 

where ip n (x) solves the non-homogeneous problem (14), 
while ^i"' 2 ^(a;) solves the problem, 

- (^ 2) )" + V{x)4? l2) - ^ l2) = X%Hn - <n& 

(25) 

The nonlinear equation (22) is just the stationary 

( 2) 

NLS equation, which has sech-solitons if sign(/i„ ; ) = 

(2) 

sign(x„ ) = sign(A„). For the focusing nonlincarity, 
a = — 1, the sech-solitons bifurcate from all band edges, 

(2) 

where [i n < 0, such that A„ < 0. It follows from 
(12) that branches of gap solitons detach from all lower 
band edges downwards the corresponding band gaps [see 
Fig. 1(a)]. For the defocusing nonlinearity, a = +1, 
the sech-solitons bifurcate from all band edges, where 
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FIG. 2: Bifurcations for the on-site and off-site gap solitons 
in a self-focusing medium (a = —1). Top: the soliton power 
Pit 1 ) = ^>{ X \ l*)dx versus fi. Solid: solitons centered at 
xo = 0, dashed: centered at xo — d/2. (a-f): spatial profiles 
of gap solitons corresponding to marked points in the upper 
plot; shading marks the minima of the potential V(x). 



(2) 

/in > 0, such that A n > 0. Therefore, branches of gap 
solitons detach from all upper band edges upwards the 
corresponding band gaps [see Fig. 1(b)]. Branches of gap 
solitons are shown on Fig. 2 for a = — 1 near the band 
edges ju, = fj,o, £t = A*3j an d M = A*4 an d on Fig. 3 for 
a = +1 near the band edges \i = fi\ and fi = fi 2 - The 
families of gap solitons have been found by solving the 
nonlinear eigenvalue problem (17) with a standard relax- 
ation technique [37]. 

The sech-solitons of the nonlinear equation (22) are 
written explicitly in the form: 



A n (X) = a n sech(n n X) 7 
where a n and k„ are found from equations, 

t-±n K-nH-n ~ u i Xn a n ZK nh l n — u 

,( 2 ) 



(26) 



(27) 



provided that sign(/A ) = sign(x« ) = sign(A„). The 
sech-type soliton envelopes (26) always have single- 
humped profile. Since A n (X) is the envelope of ip n {x), 
the resulting nonlinear bound state <$> s (x) has the oscil- 
latory structure near the band edge fj, s = \i n . 



V. BRANCHES OF GAP SOLITONS DUE TO 
SYMMETRY BREAKING 

The absence of translational invariance along the x di- 
rection, associated with the presence of the periodic po- 
tential, has an important effect on the soliton properties. 
For example, it was found that discrete solitons, bifurcat- 
ing from the first band, can be centered at (on-site) or 
in-between (off-site) potential wells. In this Section, wc 



FIG. 3: Bifurcations for the on-site and off-site gap solitons in 
a self-defocusing medium (a = +1). Notations are the same 
as on Figure 2. 



demonstrate that two branches of gap solitons, bifurcat- 
ing from all the bands, are centered at different positions 
in the periodic potential. 

The gap soliton $ s (x) near the band edge fi s = fi„ is 
represented by the perturbation series expansions (18)- 
(20), provided that the formal series converges. Parame- 
ter xo in the "slow" coordinate X = e(x — xo) determines 
the location of the bound state A n (X) with respect to the 
Bloch function ip n (x). We will show that only two values 
of xo on the period of x secure convergence of the formal 
series, in the general case. Our analysis is equivalent to 
the construction of the Melnikov function, which gives 
the distance between scparatrices in the nonlinear oscil- 
lator with a small, rapidly varying force [38, 39]. Zeros 
of the Melnikov function indicate values of xo, where the 
scparatrices intersect, so that a homoclinic orbit for the 
gap soliton exists in the nonlinear problem (17) with the 
periodic potential V(x). 

We will derive the Melnikov function [38, 39] with a 
simple but equivalent method. Derivative of the nonlin- 
ear equation (17) in x results in the following third-order 
ODE: 

+ V(x)& s - fi s <5>' s + 3cr$ 2 $' s + V'(x)<f> s = 0. (28) 

If the gap soliton $ s (x) exists, then it satisfies zero 
boundary conditions as |x| — > 00. Multiplication of (28) 
by $ s (x) and integration it over x result in the following 
constraint: 



M s (x ) 



V' (x)$ 2 s {x)dx = 0. 



(29) 



The function M s (xo) is the Melnikov function for exis- 
tence of homoclinic orbits [38, 39]. The constraint (29) is 
always satisfied if the gap soliton Q s (x) and the potential 
V(x) are symmetric with respect to the location of the 
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central peak at x = xq, such that ^ 2 s {x— x ) — & 2 (xo—x) 
and V'(x — Xq) = — V'(xq — x). More precise information 
on the constraint (29) can be obtained near the band edge 
fJ-s = fJ-m where the perturbation series expansions (18)- 
(20) are valid. The function $ e (x; X) has the power series 
expansion in e, each term of which satisfies the squared- 
pcriodic boundary conditions in x, 

<P 2 (x + d;X) = <P 2 (x;X), (30) 

and the decaying boundary conditions in X, 

lim <J> e (x;A) = 0. (31) 

We shall prove that M s (xq) is exponentially small in 
terms of e. To do so, we rewrite equation (28) for 
& e (x;X), multiply it by & e (x;X) and integrate the re- 
sulting equation over x G [0, d]. Using the periodic 
boundary condition (30), we derive the relation, 



V'(x)<S> 2 e (x;X)dx 



(32) 



Using the decaying boundary condition (31), we prove 
that 



dX V'(x)$ 2 (x;X)dx = 0. (33) 

oo JO 



As a result, the function V'(x)& 2 (x; X) is expanded in 
Fourier series in x as 

oo 

V'(x)$> 2 (x;X) = W n , m (X;e)e~, (34) 



m— — oo 



such that W n ,-m(X] e) = W n , m (X; e) and 

W n , (X;e)dX = (35) 



at any order of e. The Fourier transform of F(X) is 
defined by the standard integral: 



F(k) 



F{xy kX dx. 



(36) 



The Melnikov function M s (xq) is then expanded with the 
use of the Fourier series (34) and the Fourier transform 
(36) in the form: 



M S (X Q ) = € J2 W n ,r 



2itm 



;e e * 



(37) 



At the leading order, we have W n , m (X; 0) = A 2 l (X)wn}n, 
where u>i°L are coefficients in the Fourier series, 

oo 

V'(x)^ 2 n (x) = (38) 



The zero-order term (m — 0) in the series (37) is zero 
at any order of e, since the constraint (35) results in the 
condition: W nt o(0;e) = 0. The higher-order terms with 
larger values of |m| are exponentially smaller compared 
to the terms with smaller values of \m\ in the limit e — > 0, 
since A^k) is exponentially decaying in k. Therefore, us- 
ing exponential asymptotics, we truncate the series (37) 
by the first-order terms (to = ±1) in the limit e — > 0: 



M s (x ) = eAi cos 



27TX() 



+ ^g(w^ 1 ))+E u (39) 



where 



2tt 



and 



£ .=°(^. ( S))+°«S 



Assuming that Ai 7^ 0, we conclude from (39) that there 
are precisely two families of gap solitons bifurcating from 
two roots of the function cos(2lp-) on the period of Xq. 

We now prove that, for the squared sine potential (2), 
the values of arg(w^°{ ) are the same for all band edges as 

&rg(w { °\) = -f . It is clear from (2) that V{-x) = V(x) 
and V'\x) > for < |x| < d/2, while all Bloch wave 
squared amplitudes are symmetric, such that tp 2 (— x) = 
t^n(x). As a result, it follows from (38) that arg(u4°{) = 
— ? and the two roots of xq occur at extremal points of 
V{x): xq — and xo = f • The former (minimum) point 
corresponds to the on-site gap soliton, while the latter 
(maximum) point corresponds to the off-site gap soliton, 
in accordance with Figures 2 and 3. 

When Ax = and W nA (^;e) ^ 0, higher pow- 
ers of e are generally non-zero in the first-order terms 
(to = ±1), such that only two branches of gap solitons 
$ s (x) bifurcate in a general case. If the potential V{x) 
is special such that W n ^\ = at any order of e 

but W Ut 2 {^j e ) 7^ 0, the leading-order terms in the se- 
ries (37) become second-order (to = ±2), such that four 
branches of gap solitons $> s (x) may bifurcate from four 
roots of the function cos(^p-) on the period of xq. We 
do not know whether any special potentials V{x) may 
exist to hold the constraint W U) i (% ; e) = at any order 
of e. 



VI. LINEAR STABILITY OF GAP SOLITONS 

Stability of solitons with respect to perturbations is 
an important problem for applications. Stable states 
act as attractors, and their excitation is weakly sensi- 
tive to noise or perturbations. On the other hand, un- 
stable states tend to undergo dynamical transformations 
due to a rapid growth of initial perturbations, and this 
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behavior may be useful, for example, for switching appli- 
cations [40]. 

We study the stability of gap solitons Q s (x) by consid- 
ering the evolution of perturbed solution in the following 
form, 



ip(x,t) 



Vst [$ a (x) + (u(x) + iw(x)) e xt 
+ (u(x) + iw(xj) e Xt 



(40) 



We substitute Eq. (40) into the NLS equation (1) and 
perform its linearization with respect to the functions 
(u, w) describing small-amplitude perturbations. Then, 
we obtain coupled linear eigenmode equations where 
(it, w) is an eigenvector and A is an eigenvalue, 



C\u = —Aw, 



Cqw = Xu. 



(41) 



Here Co and C\ are Schrodinger operators with periodic 
and decaying potentials, 

dp 

Co = -^+V(x)-» s +<j$ 2 s (x), (42) 
Ci = —^+V(x)-fi s + 3a$ 2 s (x). (43) 

We are interested in eigenvalues A, which correspond 
to the spatially localized eigenvectors (u, w) in L 2 (R, C 2 ). 
If there exists an eigenvalue A with Re(A) > 0, the gap 
soliton & s (x) is spectrally unstable. On contrary, if all 
eigenvalues have Rc(A) = 0, the gap soliton is neutrally 
stable. Neutral stability can result in spectral instability 
due to resonances, embedded eigenvalues, and bifurca- 
tions of isolated eigenvalues with Rc(A) = 0. We have 
used an approach based on the Evans function for nu- 
merical calculation of the eigenvalues, the details are pre- 
sented in Appendix A. 

The stability problem (41) is written in terms of two 
Schrodinger operators Co and C\ with periodic V(x) and 
decaying a<& 2 (x) potentials. At the band edge fx = /i„, 
where $ s (x) = 0, the two Schrodinger operators coincide 
with the operator C s : 



//,, 



(44) 



For w = iu and A = the spectral bands of the 
stability problem (41) occur at Q + fj, s G Eband, i.e. 
at O G \po - Ms, Mi - Ms] U [fi 3 - fi s ,fj,2 - Ms] U .... 
For the gap soliton bifurcating from the upper band 
edge fi s = fj, n , the parameter /i s satisfies the inequality: 
fx n < n s < fJ-n+2, while for the gap soliton bifurcating 
from the lower band edge fi s = fi n , the parameter /i s 
satisfies the inequality: Mn-2 < Ms < Mn- 

We demonstrate below that an important value which 
defines many stability properties is the energy of the spec- 
tral band, which is defined by 



where the inner product (■,■)<* is defined for periodic 
Bloch functions on the period x <G [0, d\: 

(f,g)d = / f(x)g(x)dx. 



(46) 



(u,Ciu) d + (w,C w) d , 



(45) 



It is clear that h m = 2(/i m - fj, n )(ipm,ipm)d at e = 0, 
where h m refers to the m-th band edge in the spectrum 
of C s for fi s = /i„. All spectral bands of C s , which are 
lower with respect to /i s — /i„, become bands of negative 
energy for the gap soliton <& s (x), while all spectral bands 
of C s , which are upper with respect to [M s = fj, n , become 
bands of positive energy for the gap soliton & s (x)- 

The spectrum A of the stability problem (41) is double 
because of the inversion symmetry: w = —iu and A = 
— iQ. As a result, the bands of positive and negative 
energies of the operators C s and (— C s ) may overlap in 
the coupled spectrum (41) for the same values of A. 

The spectrum A of the problem (41) transforms when 
e ^ 0. A simple and stable transformation is a shift of 
spectral bands of C s and (— C s ) along the imaginary axis 
of A to the distance — fx n \. As a result, the origin 
A = becomes isolated from the spectral bands of C s 
and (—C s ) for any e ^ 0. 

Other transformations of the spectrum A are possible 
and may result in instabilities of gap solitons. These 
transformations arc considered in Sections VII and VIII. 



VII. SYMMETRY-BREAKING INSTABILITY 
OF GAP SOLITONS 

In Sec. V, we have identified two families of on-site 
and off-site gap solitons, which have different positions 
with respect to the underlying potential. In this section, 
we demonstrate that one of these soliton families is un- 
stable with respect to symmetry breaking. They tend 
to move across the potential and eventually transform 
into their stable counterparts which have a different po- 
sition. These results generalize the previously found in- 
stability of off-site discrete solitons associated with the 
first band [40]. 

More specifically, we show that the symmetry-breaking 
instability of gap solitons is defined by the sign of M' s (xq). 
If M' s (xo) > 0, then a pair of purely imaginary eigenval- 
ues A in the stability problem (41) bifurcates from A = 0, 
and these internal modes describe oscillations of the per- 
turbed soliton around the stable position x = xq. On 
the other hand, if M' s (xq) < 0, then a pair of real eigen- 
values A bifurcates in the problem (41) and these expo- 
nentially growing instability modes characterize soliton 
motion away from the unstable location x = xq. Wc 
note that these results are valid in the vicinity of gap 
edges, where the eigenvalues A are exponentially small in 
terms of the perturbation parameter e. 

Due to the symmetry of the NLS equation (1), we have 
a non-empty kernel of the operator Co for all e along the 
family of the gap soliton < I > s (x): 

Co<$>e{x-X) =0. (47) 
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On the other hand, the gap soliton ^ s (x) in the asymp- 
totic representation (19) is parameterized by xo in the 
formal power series (20) in e. As a result, the kernel of 
the operator C\ is non-empty at all power orders of e n : 



dU e = 0(e"), 



U f = 



dX ' 



(48) 



The zero eigenvalue of C\ is destroyed beyond the powers 
of e, since the gap soliton $ s (a;) is not parameterized by 
xo, values of which are fixed by roots of the Melnikov 
function (29). We show in Appendix B that the zero 
eigenvalue of C\, associated with the eigenfunction U e (x), 
shifts according to the quadratic form: 



(U e ,C 1 U e ) = ^M' s (x ), 



(49) 



where the quadratic form is defined for decaying func- 
tions on the whole line of x: 



(/,<?) 



f(x)g(x)dx. 



(50) 



According to the standard perturbation theory [41], 
the quadratic form in (49) determines the shift of the 
zero eigenvalue of C±, associated with the eigenfunction 
U e (x). When M' s (xo) > 0, the zero eigenvalue of C\ be- 
comes positive, while when M' s (xo) < 0, the zero eigen- 
value of L\ becomes negative. We show that a small 
negative eigenvalue of C\ results in a small real positive 
eigenvalue A of the stability problem (41), while a small 
positive eigenvalue C\ results in a pair of small imaginary 
eigenvalues A. 

A small eigenvalue A = A £ , corresponding to the eigen- 
function u e (x), can be found from the problem: 

Ciu € = -X^Cq 1 ^, (51) 

or equivalently, from the Rayleigh quotient: 

2 (Ue,£lU e ) 



(u £ ,£ U t 



(52) 



The quadratic form (u c ,£ u e ) exists if ($ £ , u e ) = 0, as 
follows from (47). Since ($ e ,£/ e ) = 0(e") and C x U t = 
0(e n ) at all power orders of e n , we conclude that 



u e (x) = U e (x)+E, 



(53) 



where E e is exponentially small in terms of e. We shall 
prove that 



(u e ,£ 1 u e 



4? (*«*•> + (l 



(54) 



such that (u e , C 1 u e ) > at the leading order. It follows 
from the nonlinear problem (17) that 

C X9 t (x;X) = -2e 2e — . (55) 



As a result, we have 



dX' dX 



2 £ 2 \ dX' ' 

1 /9£e r -ld$ e 

e \dX lL ° dx 



(56) 



(57) 



Solution of the inhomogeneous problem, 

r v - d ^ 

t-0 V e — — , 

ox 

exists at all power orders of e™, since the right-hand- 
side of (57) is orthogonal to $ € at all power orders of 
e™. Therefore, the quadratic form (U t ,V e ) has a regular 
power series in e, starting with the zero-order term. Since 



2e 2 \dX' ' 



1 
4? 
1 

4^2 



(* 6 ,*e) 







dX 



(X$ 2 ) dx,(58) 



and the second term is exponentially small in e, we have 
(54) at the leading order, such that the Rayleigh quotient 
(52) is given in the leading order by: 



A 



2M' a (x ) 
e 2 ($ e ,* e )' 



(59) 



It follows from (59) that a negative eigenvalue of C\ for 
M' s (xo) < results in a small positive eigenvalue A e in 
the stability problem (41). 

The exponentially small correction of the func- 
tion M s (xq) is given by the expansion (39), where 



arg ftuj^ij = — \ for the square-sine potential (2). 

Therefore, M' s (x ) > for x = and M' s (x ) < for 
Xo = |. In the former case, the gap soliton & s (x) is lo- 
cated at the minimum point of V(x) and it has a pair of 
small imaginary eigenvalues A. In the latter case, the gap 
soliton $s(x) is located at the maximum point of V{x) 
and it is unstable with a small real positive eigenvalue 
A. Fig. 4 shows unstable eigenvalues, splitting from zero 
eigenvalues, for the branches of gap solitons with xo = i. 
We note that stability of on-site and off-site solitons can 
be interchanged in more complex potentials, such as bi- 
nary supcrlattices [21]. Additionally, stability can change 
deep inside the gap, where the asymptotic analysis is not 
applicable [42]. 

Asymptotic results for NLS solitons in the lowest 
semi- infinite band gap in the focusing case (er = — 1) 
were obtained recently by Kapitula [43] in the limit 
V(x) — * 0. Branches of NLS solitons & s (x) = <&nls(^) = 



sech( A 



Xo)) in the small periodic potcn- 



dx 



dX 



tial function V(x) are defined by zeros of the function 
M s (xq), given by (29) with $ s = $NLs(a;). Stability 
of branches of NLS solitons is defined by the derivative 
M' s {xo), such that the NLS solitons bifurcating from the 
minimum points of V(x) are stable, while the NLS soli- 
tons bifurcating from the maximum points of V(x) are 
unstable. We note that the opposite conclusion is drawn 
in Ref. [43], due to an elementary sign error. 
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FIG. 4: Eigenvalues corresponding to symmetry-breaking in- 
stabilities of gap solitons centered at xo = d/2 (shown with 
dashed lines in Fig. 2) for a = — 1. 



VIII. INTERNAL MODES AND OSCILLATORY 
INSTABILITIES OF GAP SOLITONS 

Apart from symmetry-breaking instabilities analyzed 
in the previous section, we demonstrate that gap solitons 
can exhibit a different type, so-called oscillatory instabil- 
ities. Such instabilities can occur due to a resonance be- 
tween the modes corresponding to the edges of the gap in 
which soliton is localized, as was demonstrated in the case 
of a narrow gap when the coupled-mode equations are 
applicable [24]. However, the coupled- mode theory de- 
scribed only an isolated band-gap, whereas it was found 
that oscillatory instability can occur due to resonance 
between different gaps [20, 21]. Such resonances can re- 
sult in energy redistribution between the gaps and a for- 
mation of breathing structures, as was recently demon- 
strated experimentally [10]. In this section, we present 
a systematic analysis of such instabilities, and show that 
they appear when a side-band associated with the inter- 
gap resonances falls outside a band-gap. 

Oscillatory modes and instabilities are characterized 
by eigenvalues A with non-zero imaginary part of the 
stability problem (41) for f / 0. First, we show that 
new imaginary eigenvalues A with decaying eigenvectors 
(it, w) bifurcate from the band edges A = i(/i m — (i n ) of 
the same polarity as the band edge fi s = fi n . Bifurca- 
tions of internal modes occur generally at the order of 
0(e 2 ), if n„i \in- These eigenvalues are referred to as 
the internal modes of gap solitons [22, 23], and in our 
case such modes appear due to a resonance between the 
gap edges m and n. Such resonances are possible because 
a soliton induces an effective waveguide, which can sup- 
port localized modes in other gaps [44, 45]. In Fig. 5, 
we show three modes of operator Cq supported in the 
semi-infinite gap near the edge [i = /j,q by a gap soliton 
existing in the gap near the edge /i = /i2 in the case of a 
self- focusing nonlinearity (a = —1). 

Second, we show that resonance between internal 
modes of the operator C s and the bands of the inverted 
spectrum of (— C s ) occurs if the bifurcating internal mode 
of C s becomes embedded into the spectral band of (— £ s ). 
When it happens, embedded internal modes bifurcate 
generally to complex eigenvalues A, leading to oscillatory 
instabilities of the gap soliton <E> s (x). Resonant bifurca- 



-0.23 



-0.24 



-0.25 



ZL -0-26 - 



-0.27 



-0.28 



-0.29 




FIG. 5: Linear guided modes of operator Cq in the semi- 
infinite band-gap for a gap soliton shown in Fig. 2(d). 
Left: eigenvalues marked with dots (second and third ones 
are indistinguishable within the picture scale); Right: corre- 
sponding mode profiles. 



tions of complex eigenvalues A occur generally at order 
of 0(e 4 ). 

Third, we show that the internal mode of C s may occur 
near the band edge of the inverted spectrum of (— C s ). In 
this case, bifurcations of isolated, embedded, and com- 
plex eigenvalues are all possible at the order of 0(e 2 ), 
depending on the configuration of the spectral bands of 
C s and (— £ g ). 

Finally, we show that at most one internal mode can 
bifurcate from the band edge, which is closest to the zero 
eigenvalue. This bifurcation occurs generally at the order 
of 0(e 4 ). 

We emphasize that bifurcations of new eigenvalues may 
not generally occur in higher orders of e, since the band 
edges of the spectrum of L\ and Co with $ s (a0 7^ do 
not support resonances in a generic case. Bifurcations 
of new eigenvalues may occur far from the limit e = 0, 
when the spectral bands of the linearized operator get 
additional resonances at the band edges or in the interior 
points. Bifurcations of the existing eigenvalues may also 
occur far from the limit e = 0, if the existing eigenvalues 
coalesce with each other or with the spectral bands. 



A. Non-resonant bifurcations of internal modes 

Let n be the index of the band edge fi s = fi n where 
the gap soliton <f> s (%) bifurcates from. We consider a 
different band edge of the stability problem (41) with 
A = i(fi m — /i„), such that m 7^ n. We assume that 
the m-th band edge of the spectrum of C s is located 
in a band gap of the inverted spectrum of (— C s ), such 
that 2fi n ~ fi rn Sband- Using the same perturbation 
series expansions (18) and (20), we expand solutions of 
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the stability problem (41) in the perturbation series: 

u = B m (X)^ m {x) + eB' m {X)^(x) 
+e 2 u^(x,X) + 0(e 3 ), 



B m {X)i, m {x)+eB' m {X)^g\x) 
+ e 2 w%{x,X)+0{e 5 ) 



(60) 
(61) 



and 



A = i [ix m - Ms + e 2 a m + 0(e 4 )] , (62) 



(2) (2) 

where the second-order correction terms {u m , win ) solve 
the non-homogeneous system: 

£ s u$ + (fi s - v^wW = B'i U m + 2 ty&y 

+fl m B m ifj m - SaAlBjniplipm, (63) 
C s w„f + (fx s - (i m )u$ = B[' n U m + 2 (V^)') 

+fl m B m ip m - aA^Bmiplipm. (64) 

Under the constraint that (2fi s — fj, m ) ^ Sband, the 

second-order corrections («m , ) are periodic or anti- 
periodic functions of a;, when a single Fredholm condition 
is satisfied. The Fredholm condition takes the form of the 
eigenvalue problem for Q m : 



^B'; n + 2 X ^A 2 n {X)B m - n m B m = 0, 



where 



,(2) 



(65) 



(66) 



We note that Eq. (65) describes the linear modes sup- 
ported by a soliton-induced waveguide in other gaps. The 
linear problem (65) is a Schrodinger equation with the 
solvable potential (26). There is at least one isolated 
eigenvalue if sign(/x&^) = sign(xnm) = sign(fi TO ). In this 
case, the lowest eigenvalue and cigenfunction of the prob- 
lem (65) can be found explicitly as 



o _ «.2 (2) 2 



B,, 



sech Sm ( K „X), (67) 



where s m solves the quadratic equation, 



i{s m + 1) 



4// 2) v (2) 

(2) (2) ■ 
f-t-rn "X_n 



(68) 



Isolated eigenvalues J7 rn of the problem (65), when exist, 
correspond to internal modes A in the perturbation se- 
ries (62), bifurcating in the band gaps of the operators C B 
and (— jC s ) from the band edge A = i(fi m ~ Mn)- When 
sign(/x 



(2)> 



sign(xnm), the linear problem (65) does 

(2) 



not have any isolated eigenvalues. Since sign(x„m) 



(2) 

sign(x„ ) = sign(cr), we notice that all band edges 
H = Hn that support bifurcations of gap solitons in the 
nonlinear problem (17), support also bifurcations of in- 
ternal modes A in the spectrum of a selected n-th gap 
soliton. In the focusing case, a = — 1, all lower band 
edges generate internal modes A downwards the corre- 

(2) 

sponding band gaps, i.e. fi nl < and f2 m < 0. In the 
defocusing case, a = +1, all upper band edges generate 
internal modes A upwards the corresponding band edges, 

(2) 

i.e. /im > and fi m > 0. 

It is surprising that more than one internal mode A 
could be generated near the band edge A = i(/U m — n n ). In 
the case of no periodic potential V(x) = 0, perturbations 
of NLS solitons generate at most one internal mode [22, 
46]. On the other hand, perturbations of gap solitons in 
the coupled-mode equations (derived for small V(x) in 
the narrow band gap of the first resonance) may generate 
several internal modes and complex eigenvalues [24, 47] . 
In the case of finite potential V(x), the number of internal 
modes depends on the depth of the squared sech potential 
in the eigenvalue problem (65), which is determined by 

parameters of the band curvatures and /JhV and by 
the nonlinearity coefficients and Xnm- 



B. Resonant bifurcations of complex eigenvalues 

According to the general expression (40), eigenvalues 
A with non-zero imaginary part describe soliton oscilla- 
tions, which are associated with the appearance of two 
side-band spatial frequencies [i + Im(A) and /i — Im(A). 
Gap solitons are spectrally stable for small values of 
e ^ with respect to a particular resonant oscillation 
if both of the side-bands fall inside the gaps of the lin- 
ear spectrum, whereas an oscillatory instability may arise 
when one side-band appears inside a linear transmission 
band [20, 21]. This general behavior is illustrated in 
Fig. 6, where the real part of the eigenvalue is non-zero 
indicating a presence of the oscillatory instability when 
the lower side-band is inside the transmission band of 
the inverted spectrum. However, the instability is sup- 
pressed when the side-band moves inside the band gap. 
The instability shown in Fig. 6 appears due to a resonant 
coupling between a gap soliton marked "d" in Fig. 2, and 
its own fundamental guided mode in the first gap shown 
in Fig. 5. The characteristic profiles of instability modes 
are presented in Fig. 7. The top row shows the perturba- 
tion u+iw, which corresponds to higher spatial frequency 
fx + Im(A) according to Eq. (40), and we indeed see that 
this mode closely matches the guided mode profile [cf. 
Fig. 5] in agreement with the asymptotic expressions (60) 
and (61). On the other hand, the bottom row of Fig. 7 
shows the low- frequency component, which describes the 
radiation waves emitted by the soliton when (i — Im(A) 
is inside the transmission band. We note however that 
the rate of such radiation may be very small, allowing for 
the existence of long-lived oscillating, or breathing states, 
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FIG. 6: Eigenvalues corresponding to a resonance of a gap 
soliton (marked "d" in Fig. 2) with its fundamental guided 
mode in the semi-infinite band-gap. 




FIG. 7: Profiles of linear modes corresponding to a resonance 
in Fig. 6. 



such as shown in Fig. 8. Similar effects may occur due to 
a resonance with higher-order guided modes, as shown in 
Figs. 9-11. One important difference is that the associ- 
ated breathing states can have different symmetries for 
various excited modes, cf. Figs. 8 and 11. 

In mathematical terms, stability requires that all inter- 
nal modes detaching from the band edges A = i(n m — 
reside inside the band gaps of the inverted operator 
(—£,,), and the zero eigenvalue of C\ shifts to small 
imaginary eigenvalues A. When an internal mode is em- 



50 
40 
30 
20 
10 





40 -20 



20 40 



FIG. 8: Evolution of a soliton perturbed by a linear mode 
corresponding to Fig. 7. 
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FIG. 9: Eigenvalues corresponding to a resonance of a gap 
soliton (marked "d" in Fig. 2) with its higher-order guided 
mode in the semi-infinite band-gap. 



bedded into a spectral band of the inverted operator 
(— C s ), oscillatory instability of the gap soliton & s (x) 
may arise for small values of e ^ 0. Embedded imagi- 
nary eigenvalues A are known to be structurally unsta- 
ble with respect to small perturbations and, provided 
that their energy is opposite with respect to the energy 
density of the spectral band, they bifurcate into com- 
plex eigenvalues A [48]. By construction, resonance of 
internal modes of C s with spectral bands of (— L s ) is 
only possible if the internal mode, detaching from the 
band edge A = i(p rn — fi s ), has the opposite energy sig- 
nature with respect to the energy signature of the in- 
verted spectral band [i r — 2/i„ — \i m G Eband, such that 
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FIG. 10: Profiles of linear modes corresponding to a resonance 
in Fig. 9. 
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FIG. 11: Evolution of a soliton perturbed by a linear mode 
corresponding to Fig. 10. 



A = i(/x m — jtx n ) = — /i r ). Therefore, all embedded 
imaginary eigenvalues in the linearized stability problem 
(41) are expected to bifurcate to complex eigenvalues A 
in a generic case. 

We prove in Appendix C that, provided that fi r = 
2/U n - Mm € Sband, we have 



Re(A) = e 4 r m + 0(e 5 ), T m > 0, 



(69) 



where A is the eigenvalue of the bifurcating internal mode, 
given by (62). In a generic case, when T m ^ 0, the em- 
bedded imaginary eigenvalue A bifurcates to the unstable 
domain Re(A) > and leads to oscillatory instabilities of 
the gap soliton $ s (x). 



C. Marginal bifurcations of internal modes and 
complex eigenvalues 

A marginal case between non-resonant and resonant 
bifurcations occurs when internal modes detaching from 
the band edge A = i{fJ, m ~ Mn) are located in the neigh- 
borhood of the band edge A = i(/i„ — jik) of the inverted 



spectrum. We assume here that [L m , [ik, and /i s satisfy 
the resonance condition within the mismatch of order 
0(e 2 ): 



(70) 



In this marginal case, we expand the eigenvalue A and the 
eigenfunction (u, w) of the linearized stability problem 
(41) in the modified perturbation series, 



and 



where 



-u^ k (x;X)+eu^ k (x;X) 
w^l(x;X)+ew^ k (x;X) 



+ e^l(x;X) + 0(e s ) 



A = i [nm - Ms + e^rni + 0(e 4 )] , 



(71) 
(72) 

(73) 



(0) 
ink, 

(i) 

77 1 k 

(0) 
m A; 

(i) 



B m {X)^ m {x)+C k {X)i> k {x), 
B' m (X)^(x)+C' k (X)4 1 \x), 
B m (X)^ m {x)-C k (X)^ k (x), 
B' m {X)^{x)-C' k {X)^\x). 



(2) (2) 

The second-order correction terms (u mk ,w mk ) solve the 
system: 



L {2) 
l mk 

,.(2) 



(74) 



where 
77(2) 



= B'l 



+ fl mk {B m ip m - Ck^k) ~ v mk Ck^k 

- SaAl^ 2 , (B m ip m + C k ipk) , 

G (2) = B'U (V™ + 2 - C'l (i, k + 2 

+ Cl mk (B m ljj m + C k 1pk) + VmkCk^k 

- a-Alipl (B m ip m - C k 4>k) ■ 



Because of the resonance condition (70), the second-order 
corrections (u^ k , wj^L) are periodic or anti-periodic func- 
tions of x if two Fredholm conditions are satisfied. The 
two Fredholm conditions take the form of a coupled eigen- 
value problem for £l mk : 

^B^ + A 2 n (X) (<2 x ( ±B m + x { 2 k Ck) = n mk B m , 



(75) 



vfc'l + A 2 JX) (x { nl m Bm + 2j^C k ) 

— {v-mk + ^mfc) Cfe, 
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where Xnm 
defined as 



\ 2 2i is defined in (66), while Xnmk an< ^ Xnlm arc 
is 



X2) _ 



(76) 



the second case, internal modes Sl m k bifurcate as com- 
plex eigenvalues fl m k due to Fermi golden rule. Again, 
we have oscillatory instabilities of the gap soliton <& s (x) , 
emerging from all bifurcating internal modes of C s in res- 
onance with the spectral bands of (— C s ) or vice verse. 



The coupled eigenvalue problem (75)-(75) is not self- 
adjoint and therefore the eigenvalues fl m k could be 
complex-valued. 

We assume that sign(/Xm ) = sign(xnm) such that the 
first equation (75) with C% = has at least one internal 

(2) 

mode for sign(O m fc) = sign(/im '). For convenience, we 

(2) 

consider here the defocusing case a = 1, when Xnm > 

(2) 

and i±m > 0. In this case, the internal mode of the first 
equation (75) with Ck = exists for Q- m k > 0, while 
the spectral band is located for negative values of fi m fc. 
There are two particular cases, depending on whether 
[i^ > or /4 2) < 0. 

In the case fA.' < 0, the second equation (75) with 
B, m = docs not have any internal modes, while the 
spectral band is located below the value £l m k < —v m k- 
When v m k S> 1, internal modes in the component B m 
for Vt m k > are not affected by the spectral band in 
the component Ck , since the following estimate holds for 
finite Q m k and large v m k- 



(2) 

C k = -^hinA 2 n (X)B r , 

V m k 



Oh" ■ (77) 



When the value of v m k decreases and becomes negative, 
all internal modes in the component B m become embed- 
ded into the spectral band in the component Ck- The 
embedded eigenvalues Cl m k bifurcate as complex eigen- 
values Q m k due to Fermi golden rule as in [48]. 

(2) 

In the case of [i k ' > 0, the second equation (75) with 



Bm = has at least one internal mode for 



ftmfc < 0, 



> 



nik 



while the spectral band is located above the value Q 
~ v mk- When v m k ^ ~lj all internal modes in the com- 
ponent B m for Cl m k > and those in the component 
Ck for v m k + ^mk < are located in the gap between 
the two spectral bands. The internal modes in the com- 
ponent B m are not affected by the spectral band in the 
component Ck , since Ck is small according to the expan- 
sion (77). On the other hand, the internal modes in the 
component Ck arc not affected by the spectral band in 
the component B m , since the following estimate holds for 
finite (Clmk + Vmk) and large fi mfe : 



D. Internal modes near A = 

The coupled eigenvalue problem (75)-(75) is derived 
under the resonance condition (70) between two band 
edges of operators C s and (— L s ). The resonance con- 
dition (70) is always satisfied for fi m = fif. = fi n and 
Vmk = — 2A„, when the band edge A = i(/x m — (i n ) = 
of the stability problem (41) coincides with the band edge 
Ms = °f the gap soliton $ s (x) and A„ is used in (18). 
In this case, the coupled eigenvalue problem (75)-(75) 
describes the transformation of the spectrum of the prob- 
lem (41) at e 7^ 0, when a narrow spectral gap appears in 
the spectrum of the problem (41) near the origin A = 0. 

We showed in Section VII that a pair of real or purely 
imaginary eigenvalues bifurcate from A = due to the 
broken translational invariance. We will show here that 
another pair of internal mode may bifurcate inside the 
same gap from the band edges. On contrary to the former 
bifurcation, which is exponentially small in e, the latter 
bifurcation occurs generally in the order of 0(e 4 ). 

For the case /Lt m = /Ufc = M« anc l v mk = ~ 2A„, the 
system (75)-(75) can be simplified due to the obvious 



reduction: 



fj, 



(2) 



(2) = (2) 
Xnkm aVi • 



and 



(2) 



\ln ' and Xnm 



v (2) 
X n k 



(2) 
Xnmk 



Using the variables, 



B„ 



C k 



i(B m -Ck), 



A< n > = i(ft mk ~ A„), 

we transform the problem (75)-(75) to the form: 



4 n V> = -\ {r 



C {n) w 



(79) 



(80) 



(81) 



where and £ n> arc linear Schrodinger operators 

with decaying potentials: 



,(n) 



(") _ ,.(2) 



c 



,(2) 



dX 2 

d 2 

dX 2 



A„+3xi 2) ^(X), 



(82) 
(83) 



B r; 



(2) 
Xnmk 

Qmk 



A 2 n (X)C k + 



1 



n 2 , 

rnk 



(78) 



When the value v m k increases and becomes positive, the 
gap between spectral bands disappear and all internal 
modes in the components B m and Ck coalesce or become 
embedded into overlapping spectral bands. In the first 
case, internal modes VL m k bifurcate as complex eigenval- 
ues fi m due to the Hamiltonian Hopf bifurcation. In 



where sign(A r 



sign(/xl 2) ) 



(2) 

sign(x„ ). The linear 



eigenvalue problem (81) is the linearized NLS problem 
on the real line, associated to the sech-solitons (26). The 
problem has two branches of the continuous spectrum for 
A^") <G i(— oo, — |A„|) Ui(|A n |,oo), the four-dimensional 
null space A 1 -™' = and the resonance at the band edges 
AW = ±i|A„|. A small perturbation of the decaying 
potentials in the problem (81) may result in the edge bi- 
furcation of a single pair of internal modes A'™- 1 = ±iQ( n \ 
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such that fi( n ) < |A„|, provided a certain integral crite- 
rion is satisfied [23, 47]. 

It was shown [22, 46] that the discrete NLS equation 
with a small lattice step size supports bifurcations of a 
single pair of internal modes from the band edges beyond 
the linearized NLS problem (81). In order to study these 
bifurcations, we would have to extend perturbation series 
expansions (71)-(73) to the next orders and derive the 
0(e 2 ) corrections to the linearized NLS problem (81). 
This work goes beyond the scope of the present paper. 
We only note that there is at most one pair of internal 
modes bifurcating in the narrow gap near A = 0. 

IX. CONCLUSIONS 

We have presented a systematic analysis of the exis- 
tence, bifurcations, linear stability, and internal modes 
of gap solitons in the framework of the nonlinear 
Schrodinger equation with a periodic potential. This 
model or its generalizations appear in a variety of physi- 
cal applications including low-dimensional photonic crys- 
tals, arrays of coupled nonlinear optical waveguides, 
optically-induced photonic lattices, and Bose-Einstein 
condensates loaded onto an optical lattice. In the frame- 
work of this model, we have classified branches of gap 
solitons bifurcating from the band edges of the Floquet- 
Bloch spectrum, by means of the multi-scale perturba- 
tion series expansion method. We have demonstrated 
that gap solitons can appear near all lower or upper band 



edges of the spectrum for focusing or defocusing nonlin- 
carity, respectively. For the first time to our knowledge, 
we have studied the gap-soliton internal modes and sta- 
bility of gap solitons in the framework of the continuous 
model with a periodic potential. We have demonstrated 
that the gap-soliton stability is determined by the broken 
translational invariance, as well as the location of internal 
modes with respect to the spectral bands of the inverted 
spectrum. We have shown analytically and numerically 
that complex eigenvalues of the stability problem corre- 
spond to oscillatory instabilities of gap solitons. 

The analytical results presented above are rather gen- 
eral, and they are expected to be valid for different 
types of smooth arbitrary-shaped periodic potentials. Al- 
though our numerical examples have been presented for 
the specific case of the sinusoidal potential, we expect 
that our main results can be applied to other types of sim- 
ilar nonlinear models with periodically varying parame- 
ters, such as nonlinear photonic crystals and optically- 
induced photonic lattices. 
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APPENDIX A: NUMERICAL METHOD FOR CALCULATION OF EIGENVALUES 

Eigenvalues of the spectral problem (41) provide a key information about the soliton stability. However, an accurate 
numerical calculation of complex eigenvalues describing oscillatory instabilities of gap solitons is a nontrivial problem 
even in the case of a simpler system of coupled-mode equations [49, 50, 51]. The reason for numerical difficulties 
is that different components of eigenvectors have very different localization widths. For example, the modes shown 
in the bottom part of Figs. 7,10 are much broader than the soliton width, while the modes shown in the top part 
of Figs. 7,10 have comparable width. Numerical approaches used in a number of earlier studies [49, 50, 51] were 
based on the discretization of Eq. (41), however an accurate description of weakly localized modes requires the use of 
impractically wide computational windows. It was suggested that the eigenvalues can be calculated approximately, 
and then improved using a special iterative procedure [50, 51]. In our analysis, we avoid such problems by using a 
different approach based on the Evans function formalism. This method proved to be very effective for tracing soliton 
instabilities in periodic systems [52]. 

First, we reformulate the spectral problem (41) using a different set of functions U — u + iw and W = u — iw, 

--j-2+V{x)U + (T®l{x){2U + W) = (fi s +i\)U, 
d 2 W 

--j^- + V{x)W + a$ 2 s (x)(2W + U) = (ji,-i\)W. (Al) 

The advantage of this formulation for the numerical analysis is that Eqs. (Al) become uncoupled away from the 
soliton core, at |x| — > oo. In these regions, solutions of Eqs. (Al) are found in terms of Bloch functions, and they 
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form a natural basis for representation of solutions along the whole line, 

U(x) = U 1 (x)iP+ + U 2 (x)^+, 
W(x) = Wi(a:)V>r + W 2 {x)i>2, 



(A2) 



where ipi 2 ( x ) are two linearly independent Bloch functions, found as solutions of Eq. (4) with /i = fj, s ± iX and Ui. 2 
and W\ 2 are unknown parameters. By using method of variation of parameters, we set the constraints on XJ\ 2 and 
W h2 : 



dx 



dW w ( 
— = W±(x 

dx 



dx 
dx 



- U 2 (x) 



■W 2 {x 



d^4 
dx ' 

#2 



(A3) 



r/,f 



After substituting Eqs. (A2),(A3) into Eq. (Al), we obtain a set of first-order linear differential equations for the 
amplitude functions Uj(x) and Wj(x), j = 1,2, as follows 



dU. 



— f - (-l)V^(:r)(2[/ + W)^i_ 3 lV+ 



dx 
dx 



(A4) 



(-iya^ 2 s (x)(2W + U)^/V 



where the Wronskian determinants P ± = ^(dipf /dx) — ipf(dipf /dx) are x-independent [53, Sect. 1.6]. Whereas 
Eqs. (A4) are fully equivalent to the original eigenvalue problem, they are much better suited for numerical analysis 
since Ui >2 and Wi t2 only change in the region of the soliton core, where &g(x) is non-small. The key advantage is that 
the required size of the computational window is defined by the soliton width, and does not depend on the localization 
of linear modes. 

We seek spatially localized eigenmodes, which can exist when the Bloch functions 4>i 2 {x) have complex Bloch 
wave-numbers k(fx), and according to Eq. (6), one of the Bloch functions is exponentially growing whereas the other 
one is decaying. We assume, with no loss of generality, that \ipf\ — > at x — > +00. Then, a localized mode can form 
when simultaneously 



lim (U 2 ,W 2 ) = 0, lim (Ui,W!) 

x — >+oo x — 00 

In order to satisfy the limits (A5), the following determinant must vanish: 



0. 



£(X) = Dot 



(Ku( x ) u tA x ) u iJx) u^jx) \ 

Ut u {x) U+ W (x) U 2 ~ u (x) U 2 ~ w (x) 

W+ (x) W^ u {x) W.-Jx) 

\W+ U (x) W+ W (x) W 2 Jx) W 2 ~(x)J 



(A5) 



(A6) 



where four particular solutions of Eqs. (A4) are introduced according to the limiting behavior: 



lim 

■^+00 



wt 




lim 

X— S- + 00 



u l,w 



4 

2,w 



\ 




lim 

:)■ — — OC 



( u h- ) 






u 2<u 












\ w 2 ~ u J 







lim 



u 2 ~ w 
\w 2 ~ 






M 



Then, solution of Eqs. (A4) satisfying the boundary conditions (A5) is found as Uj = pt,Ut u + pt>Uj 
-,+w+ o_ n +fxr+ 1 - - ' -J- - 



and Wj 



Pw w £v>i whcrc (p£i P 



W ! Pu 1 Pli 



is an eigenvector corresponding to a zero eigenvalue of the matrix in 



I u j.u ' 

Eq. (A6). 

The coordinate x in Eq. (A6) is arbitrary, but for numerical calculations a better accuracy is achieved when it is 
chosen at the soliton center, x = xq. The function £ (A) is called Evans function, and its zeros define the location of 
eigenvalues. We approximate zeros of £ (A) by finding minima |£(A)| along the real axis and along the imaginary axis 
with a small real part, and then using a two-dimensional minimization procedure in the full complex plane. 
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APPENDIX B: DERIVATION OF (49) 



We rewrite the derivative equation (28) in the equivalent form: 

d$ e (x;X) d$ e (x;X) 
Li t: h e£i- 



dx dX 
Using the inner product (50), we reduce (Bl) to the quadratic forms: 



V'(x)$ e (x;X). 



8X' 1 dx +e \dx' x dX 



V\x)^X) d -^l d x. 



(Bl) 



(B2) 



Using the Fourier series (34), Fourier transform (36), and the expansion (37), we reduce the right-hand-side of (B2) 
to the form, 



v'(x)^x) d -^l d x = i ± 2 -^w, 



OX 2 ^ d \ ed 1 

m— — oo 



2nm 



e e d 



On the other hand, the first term in the left-hand-side of (B2) is identically zero: 

[dX Xl ~dx-) - \ 1 ~dX'~dx ' 
which is proved from the nonlinear problem (17) as follows: 



~dx 1 ~dX ~ dX 



^ + U(x)$ e - Ms $ e + ae 2 <S> 3 ) 

$'/ + V(x)<P e - n s <P e + <re 2 <S> 3 e ) = 0. 



dXdx 

The quadratic form in (B2) is therefore rewritten in the final form (49) 



\m' s {x ). (B3) 



(B4) 



(B5) 



APPENDIX C: DERIVATION OF (69) 



We extend the perturbation series expansions (60), (61), and (62) to the higher orders in powers of e. Solving the 

(2) (2) 

system (63) and (64) for the second-order correction terms (ufn , Wm ), we write the solution in the implicit form, 



u# = B>; n (X)^(x) + 2A 2 n (X)B m (X)^(x) + A 2 n (X)B m (X)4>fri 2 Hx), 
w%> = B'; n (X)^(x) + 2A 2 n (X)B m (X)^(x) - A 2 n (X)B m (X)^ (x) 



where ipm\x) is defined in (14), while iptlm (%) an d 4>^nm\x) solve the non-homogeneous problems, 



- (V4 n ™ >)" + V{x)^ - /i m Vit 2) = X^m - <"ftl> m , 
(^ 2) )" + V{x)^ + ( Mro - 2/x s )41 2 ) = -tr^Vw 



(CI) 
(C2) 



(C3) 
(C4) 



The first equation (C3) defines periodic or anti-periodic functions £ 2 '(i), since the Fredholm condition is satisfied 
by the relation (66). The second equation (C4) defines a non-periodic function 4>£m\x), since \i r G Sband, where 
/i r = 2fi s — fi m . We use the Sommerfeld radiation boundary conditions for function (j^m^ix) in the ends of the period 
x G [0, d]: 



(^nt 2) )'(d) - ik r ^(d) = 0, (4t 2) )' (0) + »M& a) (o) = 0, 



(C5) 



where fc r = fc(/Li r ) and the dispersion relation fc(/i) is defined in the Bloch functions (6). The Sommerfeld boundary 
conditions (C5) imply that the time-dependent solution ^(x, t) of the NLS equation (1) linearized at the gap soliton 
<f> s (x)e -ZAtat takes the form of outgoing radiative waves (6) directed outwards the period x G [0,d]: 



V(x,t) - <f> s (x)e 



a + {X)<j) 1 (x)e lk - x - i ' 1 - t , x^d- 
a-{X)4>2{x)e- ik - x - i ^\ x^i)+, 



(C6) 
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where a±(X) are amplitudes of the radiative waves. The Sommerfeld boundary conditions were used recently for 
embedded solitons [54]. Since the Sommerfeld boundary conditions (C5) are not symmetric, the functions ($?m (x) are 

complex- valued. The complex- valued functions 4>nm^ ( x ) result in complex- valued corrections to imaginary eigenvalues 
A in higher orders of the perturbation series expansion (62). In order to avoid lengthy analysis of the perturbation 
series equations at the third and fourth orders of e and to capture the non-zero real part of complex eigenvalues A, 
we rewrite the linearized stability problem (41) in the form of the balance equation: 

/, \\f — - \ d ( du du _ dus dw _\ . . 

(A + A) (uw — uw) = — u— —u + w— — w . (C7) 

v ' dx \ dx dx dx dx J 

Using perturbation series expansions (60), (61), and (62), we rewrite (C7) in variables x and X. The first non-zero 
term occurs at the fourth order of e and takes the form: 

-4iRe(A)^(X)^(x) = e 4 ^x;X) + Q ^ (cg) 

where the fourth-order correction term Q±(x;X) = Q < f er \x] X) + Q±\x;X) is decomposed in a periodic function 
of x and a non-periodic function of x, the latter is given by 

Qf\x; X) = «£> (fig))' - fi£> (««)' + ^ - w$ („,£>)' . (09) 

The prime in (C9) denotes the derivative in x. Integrating (C8) over the period x £ [0, d] and over the real line of X, 

(2) (2,) 

and using the explicit representation (CI) and (C2) for Um (x;X) and win (x;X), we rewrite the balance equations 
(C8) and (09) in the form: 



4 l Rc(A) n BldX^j (J i? m dx 



2e 4 ( / A\BldX 



Anl2) (l(nl2)V _ T(„Z2) ( Anl2)\' 



+ 0(e 5 ). (CIO) 



Using the Sommerfeld boundary conditions (05), we finally derive the expansion (69), where 

The formula (Cll) is referred to as the Fermi golden rule of radiative decay of embedded eigenvalues [48, 54] . 
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